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Abstract 

Statistical methodology for the design and analysis of clinical Phase II dose response studies, 

with related software implementation, are well developed for the case of a normally distributed, 

homoscedastic response considered for a single timepoint in parallel group study designs. In 

practice, however, binary, count, or time-to-event endpoints are often used, typically measured 

repeatedly over time and sometimes in more complex settings like crossover study designs. In 

this paper we develop an overarching methodology to perform efficient multiple comparisons 

and modeling for dose finding, under uncertainty about the dose-response shape, using general 

parametric models. The framework described here is quite general and covers dose finding using 

generalized non-linear models, linear and non-linear mixed effects models. Cox proportional 

hazards (PH) models, etc. In addition to the core framework, we also develop a general 

purpose methodology to fit dose response data in a computationally and statistically efficient 

way. Several examples, using a variety of different statistical models, illustrate the breadth of 

applicability of the results. For the analyses we developed the R add-on package DoseFinding, 

which provides a convenient interface to the general approach adopted here. 
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1 Introduction 



Finding the right dose is a critical step in pharmaceutical drug development. The problem of 
selecting the right dose or dose range occurs at almost every stage throughout the process of 
developing a new drug, such as in microarray studies [25], in- vitro toxicological assays [8], animal 
carcinogenicity studies [H], Phase I studies to estimate the maximum tolerated dose [26], Phase 
II studies covering dose ranging and dose-exposure-response modeling [20], Phase III studies for 
confirmatory dose selection, and post-marketing studies to further explore dose response in specific 
subgroups defined by region, age, disease severity and other covariates [311 |32l [22] . In recent years, 
considerable effort has been spent on improving the efficiency of dose finding throughout drug 
development [371 [121 [21] . Despite these efforts, however, improper dose selection for confirmatory 
studies, due to lack of sufficient dose response knowledge for both efficacy and safety at the end of 
Phase II, remains a key driver of the ongoing pipeline problem experienced by the pharmaceutical 
industry [5l[28]. 

Statistical analysis methods for late development dose finding studies can be roughly categorized 
into modelling approaches to characterize the dose response relationship [301 ISSj and multiple test 
procedures for dose response signal detection [M] or confirmatory dose selection [351 IE!- Hybrid 
approaches combine aspects of multiple testing with modeling techniques to overcome the short- 
comings of either approach [381 [10]. More recently, considerable research has focused on extending 
these methods to response-adaptive designs that offer efficient ways to learn about the dose response 
through repeated looks at the data during an ongoing study [TOl [13, [231 [1] • 

Most statistical methodology for dose response analysis has been introduced in the context of a 
normally distributed, homoscedastic endpoint, with a parallel group design, in which each patient 
receives only one treatment. In practice, however, one often faces more complex study designs 
(e.g., cross-over designs), where a heteroscedastic or non- normally distributed endpoint is measured 
(e.g., binary, count and sometimes time-to-event data). One approach is to extend the existing 
methodology using generalized non-linear models or generalized non-linear mixed effects models. 
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However, these extensions are typically specific to each new situation. In addition general purpose 
software for these types of models is not available and a case-by-case implementation requires a 
major coding effort. In this paper, we describe an overarching hybrid approach, combing multiple 
comparisons and modeling, to the analysis of dose response data for general parametric models and 
general study designs, that allows for a straightforward computer implementation. 



Ti'ial design stage 

• Deteimioatiou of the endpoint, study population, etc. for which the dose response 
rehitionship sliould be establislied 

• Pre-specificatioii of candidate set of dose-response models 

• Sample size calculation, based on either power considerations to establish a dose 
response signal or the precision in dose response / target dose estimation 



Data collection 



Trial analysis stage 

• "MCP"step 

o Assessment of a dose response signal using a suitable test for trend 
o Model selection (or model a^ eraging) out of the set of statistically significant 
dose response models 

• "Mod" step 

o Dose response and target dose estimation based on the selected model 



Figure 1: Overview of MCPMod approach 

More specifically, we extend the MCPMod approach [10], which was originally introduced for nor- 
mal, homoscedastic, independent data. This approach provides the flexibility of modeling for dose 
response and target dose estimation, while accounting for model uncertainty through the use of 
multiple comparison and model selection/averaging procedures. The approach can be described in 
two main steps (Fig. [1]). At the trial design stage the clinical team needs to decide on the core 
aspects of the trial design, as in any other trial. Specific for MCPMod is that a candidate set of 
plausible dose response models is pre-specified at this stage, based on available pharmacokinetic 
data/dose response information from similar compounds, etc. This gives rise to a set of optimal 
contrasts used to test for the presence of a dose response signal consistent with the corresponding 
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candidate models. In case of large model uncertainty, this candidate set should be chosen to cover 
a sufficiently diverse set of plausible dose response shapes 

The trial analysis stage consists of two main steps: The MCP and the Mod steps. The MCP 
step focuses on establishing evidence of a drug effect across the doses, i.e., detecting a statistically 
significant dose response signal for the clinical endpoint and patient population investigated in the 
study. This step will typically be performed using an efficient test for trend, adjusting for the fact 
that multiple candidate dose response models are being considered. If a statistically significant 
dose response signal is established, one proceeds with determining a reference set of significant dose 
response models by discarding the non-significant models from the initial candidate set. Out of this 
reference set, a best model is selected for dose estimation in the last stage of the procedure, using 
existing non- linear regression methods [2]. The selected dose- response model is then employed to 
estimate target doses using inverse regression techniques and possibly incorporating information on 
clinically relevant effects. The precision of the estimated doses can be assessed using, for example, 
bootstrap methods. 

The original MCPMod method proposed by (TU] was intended to be used with parallel group designs 
with a normally distributed, homocedastic response. Although that covers a good range of dose 
finding designs utilized in practice, the restrictions of the original method create serious practical 
limitations to its wider use in drug development. For example, binary, count, and time-to-event 
endpoints, though frequently used in many disease indications, are not covered by the original 
MCPMod formulation. Likewise, longitudinal patient data, like in crossover studies, with its im- 
plicit within-patient correlation, cannot be properly analyzed with the original formulation of the 
MCPMod methodology. In what follows, we extend the MCPMod methodology, to perform dose 
response modeling and testing in the context of general parametric models and for general study 
designs, in a similar way as [19] extended certain simultaneous inference procedures. Note that, 
even though we focus on extending the MCPMod approach, the results of this paper remain valid, 
in particular, if only a multiple comparison or a modeling approach is to be applied. 

We introduce the proposed extension in Section |2l In Section [3] we describe an efficient approach 
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for dose response model fitting, which is evaluated in terms of asymptotic considerations and sim- 
ulations. In Section |4] we use several case studies to illustrate the proposed methodology and its 
implementation with the R add-on package DoseFinding [6]. In summary the method is illustrated 
for binary, overdispersed count, time-to-event data (based on the Cox PH model) and longitudinal 
data with patient specific random effects. 

2 Generalized MCPMod 

This section describes an extension of the original MCPMod approach that can be used whenever 
the response variable can be described by a parametric model in which one of the parameters 
captures the dose response relationship. We show how the basic ideas and concepts of the original 
MCPMod can be extended to this setting. 

2.1 Basic concepts, notation and assumptions 

In the original description of MCPMod, the expected value of the response is utilized as the parame- 
ter capturing the dose response relationship. The key idea of the extended version of MCPMod is to 
decouple the dose response model from the expected response, focusing, instead, on a more general 
characterization of dose response via some appropriate parameter in the probability distribution of 
the response. To be more concrete, let y denote the response vector of an experimental unit in the 
trial (e.g., a patient) which has been assigned a dose x. The following results can easily be extended 
to the case of multiple doses if needed. We assume that y follows the distribution function given 
by 

y F {z,r],fi{x)) , (1) 

where fi{x) denotes the dose response parameter, r/ the vector of nuisance parameters, and z 
the vector of possible covariates. The key features of extended version of MCPMod can then be 
formulated with respect to /x(a;), such as: 
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• accounting for uncertainty in the dose response model via a set of candidate dose response 
models, 

• testing of dose response signal via contrasts based on dose response shapes, 

• model selection via information criteria, or model averaging to combine different models, and 

• dose response estimation and dose selection via modeling. 

Because all dose response information is assumed to be captured by the interpretability of 

this parameter is critical for communicating with clinical teams, choosing candidate dose response 
shapes, specifying clinically relevant effects for target dose estimation, etc. To better illustrate this 
point, consider a time-to-event endpoint that is assumed to follow a WeibuU distribution. The 
Weibull distribution is typically parameterized by a scale parameter A and shape parameter a, 
neither of which has an easy clinical interpretation. For the purpose of MCPMod modelling the 
model could be reparameterized in terms of the median time to event /i = [log(2)]^^" /A and a, and 
then one would use as an interpretable dose response parameter. 

All dose response models of interest in this paper can be expressed as 



where denotes the so-called standardized model function and 6 its parameter vector. For 
example, for f^{x,d^) = x one obtains the linear model f{x,6) = 60 + 9ix and for f^{x,6^) = 
x/{x + 6^) the Emax dose response model f{x, 6) = 6q + 6ix/{x + 6^)] see [IOl[^ for more examples. 
Dose response models of the form ([2]) are specified as candidate models for Covariates may 

also be included in (|2]), at the model fitting stage, but we leave them out for now to keep the 
notation simple. 

Assume that K distinct doses xi, . . . ,xk are utilized in a trial, with xi denoting placebo. Assume 
further that M candidate models /i, . . . , J'm have been chosen to capture the model uncertainty 
about We define the dose response parameter vector associated to candidate model m as 

fJ-m = {f^m,l, • • • , , whcrC firn,i = fm{xu 6) , i = 1, . . . , K, ITl = 1, . . . , M . 




(2) 
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For the purpose of obtaining estimates of the dose response, we initially consider an analysis-of- 
variance (ANOVA) parametrization for the dose response parameter = /Xj, i = 1,...,K. 

That is, we allow a separate parameter /ij to represent the dose response at each dose level. Let /x 
denote the vector of estimated dose response parameters under the ANOVA parametrization, ob- 
tained using the appropriate estimation method for the general parametric model ([1]) via maximum 
likelihood (ML), generalized estimating equations (GEE), partial likelihood, etc. Note that these 
type of ANOVA estimates are easily available from standard statistical software packages. The key 
assumption underpinning the extended version of MCPMod is that p, has an approximate distri- 
bution A^(/x, S) , where S denotes the variance-covariance of /x. This assumption can be shown 
to hold for most parametric estimation problems, such as, generalized linear models, parametric 
time-to-event models, mixed-effects models, GEE, etc. Note that S is a function of n and converges 
to as n — 7- oo. Furthermore, S may or may not depend on x. For example, dependence on x 
arises if unequal variances for different dose levels Xi are considered. The estimation of is done in 
a separate second stage based on p, and an estimate S of its covariance matrix. Section [3] explains 
this in detail. 

2.2 Implementation of the MCP step 

The MCP step consists of specifying a set of candidate models for the dose response relationship 
/i(x). To that end, one needs to specify families of candidate models {e.g., linear, Emax, logistic, 
quadratic). In addition, to derive optimal model contrasts, one needs to determine guesstimates 
for the non-linear parameters 6^, such as the ED50 parameter for the Emax model. Note that the 
shape of the dose response model function is determined by the parameter 6^, which is why only 
guesstimates for this parameter are needed to derive optimal model contrasts. As mentioned earlier, 
the clinical interpretability of fi{x) is critical for this step. Further details on and strategies for the 
specification of candidate models and corresponding guesstimates are given in [29] . 

Given these guesstimates, each candidate model shape determines an optimal contrast for a trend 
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test to evaluate the associated dose response model signal, such as a linear trend or a trend based 
on an Emax model with ED50 = 2. The optimal contrasts are applied to the previously described 
ANOVA estimates p,, with the associated asymptotic distribution used for implementing the cor- 
responding tests [i.e., critical values and p- values). It can be shown that the (optimal) contrast 
for testing the hypothesis of a flat dose-response profile with maximal power for a single candidate 
model shape /x^ is given by 

oc 5- (mII - , (3) 

where fj!^^ = (/^(xi, 0^), . . . , f^{xK, 0^))' and 6^ is the parameter for which guesstimates are re- 
quired, see [TU] and Appendix A. For convenience, we normalize the contrast coefficients such that 

1 1 f^opt II 1 • 

The implementation of contrast tests for the candidate models is done similarly to the original 
MCPMod approach. Let c°^*, . . . , c^^* represent the optimal contrasts corresponding to the set 
of candidate models and = [c°^* ■ ■ ■ c^^*] the associated optimal contrast matrix. The con- 
trast estimates are then given by (C°^*)'/2, being asymptotically normally distributed with mean 
(C^*) V and covariance matrix (C°^*) 5C°^*. It follows that the asymptotic z-test statistic for 
the m*^ candidate model hypotheses Hq : (c^*)'/x = vs. Hq : (c^*)'/x > is given by 
Zm = (c^*)'/i/ (C°^*) SC°P^ , with denoting the m*'' diagonal element of the ma- 

trix A. The test statistic used for establishing an overall dose response signal is the maximum 
Z(^M) = niaXm Zm of the individual model test statistics. Critical values for tests with (asymptoti- 
cally) exact level a can be derived from the joint distribution of z = {zi, . . . , z^), which is easily 
obtained from the joint distribution of the contrast estimates given previously, and using 

Piz^M) >q) = l- Piz^M) <q) = l-P{z< ql). (4) 

Multiplicity adjusted p-values for the individual model contrast tests can be derived similarly. The 
mvtnorm package in R includes functions to calculate quantiles and probabilities for the underlying 
multivariate normal distributions [T7] . 

If the optimal contrasts and the critical values also depend on one needs guesstimates for nuisance 
parameters appearing in the covariance matrix at the planning stage, as well. This is a difference 
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compared to the normal homoscedastic setting without covariates, where S is proportional to a 
diagonal matrix with the the reciprocal of the group sample sizes on the diagonal. Once data 
becomes available, the z-statistics for the model contrast tests are calculated and their maximum 
used for the dose response test. At this stage, one can obtain the estimated S matrix from the 
observed data, and use this to recalculate optimal contrasts and the critical value for the test. Note 
that the guesstimates for the parameters 9^ are not recalculated based on the observed data, as 
this would lead to a serious Type I error rate inflation. For the purpose of the multiple contrast 
test, the guesstimates pre-specified at the planning stage for 9^ should be used. 



2.3 Implementation of the Mod step 

Once a dose response signal is established, one proceeds to the Mod step, fitting the dose response 
profile and estimating target doses based on the models identified in the MCP step. There are 
many ways to fit the dose response models to the observed data, including approaches based on 
maximizing the likelihood (ML) or the restricted likelihood. However, a direct ML approach requires 
the derivation of the likelihood in every specific case, with a considerable amount of model-specific 
coding involved. We therefore suggest an alternative two-stage approach to dose response model 
fitting that utilizes generalized least squares. This approach is described in more detail in Section El 
It relies on asymptotic results, but has the appeal of being of general purpose application, as it 
depends only on /i and S. In addition, as shown later in the simulation study, its finite and large 
sample properties are similar to those of the approaches based on the full likelihood. 

Model selection can be based on the maximum z-statistics, or using information criteria, such as 
the AIC or the BIG. Estimates of the latter under the model fitting approach are discussed in the 
next section. Estimation of target doses is done based on the selected fitted model for the dose 
response parameter [7]. 

Alternatively, model averaging approaches can be used to avoid the need to select a single model. 
The individual AIC and BIG values for the candidate models with significant contrast test statistics 
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determine the model averaging weights [TT]. This apphes both to dose response and target dose 
estimation. 

The generahzation of MCPMod described in this section has focused on testing and estimation 
associated with the full dose response profile, that is, including the response at placebo and the 
entire dose range utilized in the study. In practice, there are cases in which inference might focus 
on the placebo-adjusted dose response (or more generally a control-adjusted response), that is, 
the dose response with the placebo or control effect subtracted fc{x,0) = f{x,0) — f{O,0). This 
could become relevant, for example, if covariates are added to the placebo response 6'o in (|2]). In 
the context of time-to-event data, the focus on placebo-adjusted dose response will occur naturally 
when modeling the hazard ratio as a function of dose. As shown in Appendix |Bl all results presented 
in this section, including the derivation of optimal contrasts, as well as the model fitting results 
described in the next section, apply equally in the context of placebo-adjusted dose response. 



3 Non-linear dose response model fitting using a two-stage, 
generalized least squares approach 

In this section we describe an efficient methodology for fitting nonlinear dose response models that 
can be used for the Mod step of the MCPMod methodology. The fitting is done in two stages: 
First, the ANOVA estimates and S introduced in Section |2] are obtained. Second, the parameter 
d is estimated by fitting the dose response model to the ANOVA estimates from the first step 
using a generalized least squares estimation criterion. In Section 13. 1[ we establish consistency and 
asymptotic normality of this estimate. In Section 13.21 we assess the accuracy of the asymptotic 
results via a simulation study. 
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3.1 Asymptotic results 



Assume that we have dose response estimates p, = (/ii, . . . jJIr)' obtained from an ANOVA-type 
parameterization of /x which allows for unrelated mean responses at each of the K dose levels; see 
Section [2Tn These estimates are assumed to be asymptotically multivariate normal distributed with 
a covariance matrix consistently estimated by S. The estimates p. and S are easily available from 
standard statistical packages, see Section H] for examples using R. Next, we fit the non-linear dose 
response model f{x,6) to the estimates p, by minimizing the generalized least squares criterion 

^{6) = {n- fix, e))'A^{p - fix, 6)) (5) 

with respect to 9 to obtain the estimates 6. In Equation (|5]) we have f{x, 6) = {f{xi, 6), . . . , f{xK, 0))' 

p p 

and An denotes a symmetric positive definite matrix. We assume that A„ — A, where — )■ denotes 
convergence in probability. In practice we will always use An = S , but this would unnecessarily 
restrict the discussion at this stage. 



Let 6q denote the true value of the parameter 6. In Appendix Owe show that, under mild regularity 
conditions, is a consi 
multivariate normality 



conditions, is a consistent estimator of O^, i.e., 6 — )■ 9q. Furthermore, we have the asymptotic 



^(9 - 0o) A N{0, B{eo)'M{e^)B{e^)), (6) 

where M(6>) = anF{e)'ASAF{e) and B(6>) = (F(6>)'AF(0))-\ F(6>) denotes the d x A; matrix 
of partial derivatives '^'^^^^f'^ {i = l,...,/c, h = a„ a non- decreasing sequence of values 

increasing to infinity as the sample size n goes to infinity. 

Selecting An = S would be the best choice, if S were known. Provided that anS — )■ 5] as n — )■ oo, 
the previous formulas in this case simplify to ^{6 - Oo) 4 N{0, (F(0o)S"^F(0o)')"^) with 
minimizing 

= (m - fix, e)ys-\n - fix, 6)) (7) 

with respect to 0. Because S is not known, we will typically use a consistent estimate S of it in ([7]) . 
If the assumptions about the covariance matrix are wrong in the sense that a„S does not converge 
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to S, then will still be consistent and have the asymptotic normal distribution given in ([6]). In 
this regard, the suggested estimator is similar to Huber's robust estimator ([21]; Section 4.6 in [13]). 
but this aspect will not be further utilized. 

Note that this two-stage, generalized least squares estimate is quite similar to the ML estimate: For 
normal homoscedastic data both approaches lead to exactly the same estimate, while, for example, 
in generalized linear model settings the two estimators have the same large-sample variance; see 
Chapter 4.3 in [H]. 

The computational advantage of using this two-stage approach is that the target function in ([7]) 
that is optimized numerically is low-dimensional: The dimension is equal to the number of different 
dose levels and the target function can thus be evaluated quite efficiently, while the target function 
in a full likelihood approach depends on the complete data set. This difference in speed becomes 
relevant in clinical trial settings, as often extensive clinical trial simulations are used to evaluate 
proposed study designs. Another advantage of the proposed method is its broad applicability to 
general parametric models. 

Model selection criteria are generally defined as — 21og(L) -|- dim(0)r, where L denotes the likeli- 
hood function evaluated at the maximum likelihood estimate and r a penalty for the number of 
parameters, which depends on the model selection criterion. One approach to model selection is 
hence to use ^{0) + dim(0)r to compare different dose response models fitted based on the same 
p, and S. This criterion is motivated by the fact that, for normally distributed homoscedastic data 
without covariates, these two approaches are equivalent in terms of selecting the same model: The 
likelihood function can be split into the sum of the deviation between the observed data and p, and 
the deviation between Jl and f{x, 6). The deviation of the individual data and p, is identical across 
the different dose response models, so that the criterion only varies with the deviation between p 
and f{x,6), which, in case of normal data, is equal to ^{0). In situations beyond the normal 
case both approaches might lead to slightly different results, however as ([7]) is roughly proportional 
to — 2-log-likelihood of the ML estimate of 0, when discarding the contribution of the first stage 
ANOVA-type fit (which is equal for all dose response models considered) typically both approaches 
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will lead to very similar results. Subsequently r = 2 will be used and we refer to this criterion as 
gAIC. 

In the next subsection we investigate the properties of the proposed asymptotic approximations via 
simulations before illustrating it with several applications in Section HI 

3.2 Simulations 

In this section we evaluate the asymptotic performance of the approximations provided in Section 
13.11 for fitting a single nonlinear dose response model. More specifically, we compare the proposed 
methods with more traditional maximum likelihood estimation by evaluating the dose response 
estimation accuracy using simulations. In addition we assess the coverage probability of the resulting 
confidence intervals for the model parameter 0. 

3.2.1 Design of simulation study 

Throughout the simulations we assume five active dose levels 0.05, 0.2, 0.5, 0.8, 1 plus placebo. We 
investigate equal group sample sizes of 15, 30, 50, 100, 300 and 1000 patients for each dose. The 
lower range of the investigated sample sizes is realistic for typical Phase II dose response studies, 
while the sample sizes of 300 and 1000 are included to assess the asymptotic behavior. 

We investigate three types of data: binary data, overdispersed count data using a negative binomial 
model and time-to-event data using a Cox PH model for estimation. Regarding dose response 
models, we will utilize an Emax model fi{x,6) = + 0ix/{02 + x), a quadratic model fi{x,6) = 
6o + Oix + 6'2X^ and an exponential model fi{x,6) = Oq + 6i{exp{x/62) — 1). In the simulations, 
we set 6*2 = 0.05 for the Emax model, 62 = 0.2 for the exponential model and 61/62 = —5/8 for 
the quadratic model; see Figure [2] for the underlying model shapes. The remaining parameters 6q 
and 61 are chosen such that the power for testing the dose with the maximum treatment effect 
against placebo at the 5% one-sided significance level is 80% for 30 patients per group. This 
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Figure 2: Dose Response Models used for simulation 

ensures a realistic range of sample sizes (in terms of the signal to noise ratio) is investigated in the 
simulations. 



Data type 


Quadratic 


Emax 


Exponential 


binary 
count 

time-to-event 
normal 


(-1.734,4.335,-2.7094) 
(2,-2,1.25) 
(0,-1.8876,1.1797) 
(0,2.61,1.633) 


(-1.734,1.8207,0.05) 
(2,-0.84,0.05) 
(0,-0.7928,0.05) 
(0,1.097,0.05) 


(-1.734,0.01176,0.2) 
(2,-0.005427,0.2) 
(0,-0.005122,0.2) 
(0,0.007089,0.2) 



Table 1: Dose Response Model Parameters {0o,9i,62) used for simulation 



Table [T] summarizes the three dose response model specifications for each data type. For binary 
data. Table [1] gives the the mean on the logit scale. For count data, the logarithm of the mean 
is as specified in Table [H and the overdispersion parameter is 1. For time-to-event data, we use 
an exponential distribution for data generation with the log-means specified in Table [1] and where 
observations larger than 10 are censored. The mean in the placebo group is 1, so that the log-mean 
is equal to the difference in log-hazard rates. The Cox PH model is formulated relative to the control 
group, so that, in this case, the placebo parameter is set to when estimating the dose response 
model. In addition, we include normally distributed data with a = 1 as a benchmark comparison, 
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since in this case the two-stage, generahzed least squares (GLS) and ML estimates coincide. 

We use the two-stage approach from Section [2] to obtain (i) an ANOVA-type model fit to the data 
by using either a generalized linear model (binary and count data) or a Cox PH model (time-to- 
event data) with "dose" treated as a factor and (ii) a dose response model fit to the resulting dose 
response estimates obtained via generalized least squares ([7]), together with the asymptotic results 
from Section 13.11 In the simulations, we compare this approach to nonlinear maximum likelihood 
(binary and count data) and maximum partial likelihood (time-to-event data) estimation using the 
same link functions as above. For the model fitting step, we assume lower and upper bounds for the 
62 parameter of 0.001 and 5 for the Emax and 0.05 and 5 for the exponential dose response model. 

In addition, we assess the coverage probability for three different methods of constructing confidence 
intervals for the dose response model parameter 0. First, we use the generalized least squares (JTj) 
together with the asymptotic results from Section 13.11 (denoted below as GLS). Second, we use 
parametric bootstrap confidence intervals by sampling from the multivariate normal distribution 
underlying the first stage ANOVA-type estimates and then fitting the nonlinear model to each 
of these samples using the GLS criterion from ([7]). The bootstrap confidence intervals are then 
calculated by taking the 5% and 95% quantiles of the observed sample. For each simulation we 
used 500 bootstrap samples (denoted below as GLS-B). Finally, we use the maximum likelihood 
fits and calculate confidence intervals based on the inverse of the Hessian matrix and the usual 
asymptotic normality assumptions (denoted below as ML). 

3.2.2 Results of simulation study 

Simulations were run with 2000 replications, using the DoseFinding package version 0.9-1. To 
illustrate the performance of the GLS and ML methods with regard to dose response estima- 
tion, we calculated the root mean squared estimation error averaged over the available doses 
J\Y.Ui.x^O)- f{x, e)y, where V = {0, 0.05, 0.2, 0.6, 0.8, 1}. Figures [9l [lO] and [ID in the Ap- 
pendix display the results. It is evident from these plots that both approaches can hardly be distin- 



15 



guished in terms of the mean squared error, indicating that, in terms of dose response estimation, 
both methods perform almost identically, even for small sample sizes. 

Next, we assess the coverage probability for the three different methods described at the end of 
Section I3.2.1I Figure |3] displays the results only for the count data case, because the results for the 
binary and time-to-event data are nearly indentical. We observe that the asymptotic confidence 
intervals for the GLS and ML methods perform very similarly for all three models under investigation 
(Emax, exponential, quadratic). Both methods achieve the nominal 90% coverage probability fairly 
well for the quadratic model, even for small sample sizes. For the Emax model, the nominal coverage 
probability is achieved at roughly 50-100 patients per group, whereas for the exponential model the 
coverage probability is achieved only at very large sample sizes (the poor performance of standard 
asymptotic confidence intervals for nonlinear regression models even for moderate sample sizes is 
well-known, see for example, [27]). The reason for the better performance under the quadratic 
model is the fact that it is linear in the model parameters. One reason for why the confidence 
intervals perform worse for the exponential model than for the Emax model might be that the dose 
design used in the simulations allows an easier identification of the model parameters under the 
Emax model, because there are more dose levels in the lower part of the dose range than the upper 
part. 

In contrast, the parametric bootstrap approach GLS-B achieves the 90% nominal coverage prob- 
ability fairly well for all three dose response models, even at sample sizes as small as 30 patients 
per group. The GLS-B method thus performs always at least as well as the GLS and ML methods, 
and the general recommendation is to use this in case of small sample sizes. The approach is com- 
putationally more expensive, as it requires repeated fitting of the dose response models, but the 
bootstrap-based on the GLS two-stage fitting is computionally much more efficient than a traditional 
bootstrap approach based on ML: The GLS two-stage approach only depends on the ANOVA type 
estimates and not the individual observations, which makes evaluation of ^{0) computationally 
much cheaper than evaluation of the full likelihood function. 

In summary, we conclude that both the GLS and ML methods perform similarly under the different 



16 



thetaO thetal 



theta2 



30 50 100 300 1000 



1.00 
0.95 
0.90 
0.85 
0.80 



1.00 - 

0.95 

0.90 

0.85 

0.80 



GLS-B 


GLS-B 


GLS-B 


emax 


exponential 


quadratic 












































1. . 




P— si 






























^-1 


K 




r 


















^ ' 1 






r- 












1 




















































































ML 


ML 


ML 








emax 










exponential 










q 


jadra 


ic 
























' ■ ■ ■ ■ 


* . 


1. 






















■ - H 


— — 































































> 


■ 








1^ 






























r-^ — 


■ 





























































GLS 


GLS 


GLS 








emax 








exf 


)onen 


tial 










q 


jadra 


ic 




















— 1 




























- ^ 


— 






































































































































' — -A 


1 ^ 


• — -« 


1 





















1.00 
0.95 
0.90 
0.85 
0.80 



15 30 50 100 300 1000 



30 50 100 300 1000 



;ure 3: Empirical coverage probability (based on 2000 simulations), of 90% confidence intervals. 



17 



dose response shapes, sample sizes, and data types investigated in the simulation study. This 
conclusion holds both for the coverage probabilities, as well as for the average estimation error. As 
mentioned previously, however, the GLS method is very general and computationally more efficient, 
which facilitates the usage of computationally expensive techniques such as the bootstrap approach. 

4 Applications 

4.1 Longitudinal modeling of neurodegenerative disease 

This example refers to a Phase 2 clinical study of a new drug for a neurodegenerative disease. The 
state of the disease is measured through a functional scale, with smaller values corresponding to 
more severe neurodeterioration. The goal of the drug is to reduce the rate of disease progression, 
which is measured by the linear slope of the functional scale over time. 

The trial design includes placebo and four doses: 1, 3, 10, and 30 mg, with balanced allocation of 
50 patients per arm. Patients are followed up for one year, with measurements of the functional 
scale being taken at baseline and every three months thereafter. The study goals are to (i) test the 
dose-response signal, (ii) estimate the dose-response and (iii) select a dose to be brought into the 
confirmatory stage of the development program. 

The functional scale response is assumed to be normally distributed and, based on historical data, 
it is believed that the longitudinal progression of the functional scale over the one year of follow 
up can be modeled by a simple linear trend. We use this example to illustrate the application of 
MCPMod in the context of mixed-effects models. 

We consider a mixed-effects model representation for the functional scale measurement yij on patient 
i at time tij-. 

Vij = (/3o + boi) + {n{xi) + bii) tij+eij, [boi, fen]' (0, A) and e-ij^N (O, cr^) , all stoch. independent. 
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The DR parameter in this case is the hnear slope of disease progression fi{x). If fi{x) is represented by 
a hnear function of dose x, the model in ([8]) is a linear mixed-effects (LME) model, else it becomes a 
nonlinear mixed-effects (NLME) model. In particular, under the ANOVA parametrization discussed 
in Section [2TT1 fi{xi) = fii, ([8]) is an LME model with different slopes for each dose. 

The research interest in this study focuses on the treatment effect on the linear progression slope. 
At t = 1 year this is numerically equal to the average change from baseline, and thus easily 
interpret able. At the planning stage of the trial, the following assumptions were agreed with the 
clinical team for the purpose of design: 

• Natural disease progression slope = -5 points per year. 

• Placebo effect = {i.e., no change in natural progression). 

• Maximum improvement over placebo within dose range = 2 points increase in slope over 
placebo. 

• Target (clinically meaningful) effect = 1.4 points increase in slope over placebo. 

Guesstimates for the variance-covariance parameters were obtained from historical data: var (6oj) = 
100, var (bii) = 9 corr (6oi, bu) = —0.5, and var (ejj) = 9. Under these assumptions, it is easy to see 
that the covariance matrix of the ANOVA-type estimate p, of the slopes /x = (nimg, l^smg, l^iomg, l^somg)' 
is compound-symmetric. With these concrete guesstimates, the diagonal element is 0.1451 and the 
off-diagonal element 0.0092. 

From discussions with the clinical team, the four candidate models displayed in Figure H] were 
identified: 

• Emax model with 90% of the maximum effect at 10 mg, corresponding to an ED50 = 1.11 

• Quadratic model with maximum effect at 23 mg, corresponding to standardized model pa- 
rameter 6 = —0.022 
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Figure 4: Candidate models for neurodenerative disease example. 

• Exponential model with 30% of the maximum effect occurring at 20 mg, corresponding to a 
standardized model parameter 6 = 8.867 

• Linear model 

For confidentiality reasons, the data from the actual trial cannot be used here, so we utilize a 
simulated dataset with characteristics similar to the original data with an Emax DR profile imposed 
on the linear slopes fi{x). Figure E] shows the simulated data per dose, which is available in the 
DoseFinding package, in the neurodeg data set. 

In what follows we illustrate the individual steps of MCPMod along with its implementation in 
DoseFinding package (version 0.9-6). 

The p, vector is estimated via an LME fit of data, which can be done, for example, using the Ime 
function in the nlme R package, as illustrated below. 

data (neurodeg) 
head (neurodeg, n=3) 
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Figure 5: Simulated data for neurodegenerative disease example. Gray lines correspond to individual 
patient profiles, black line to loess smoother. 

> resp id dose time 

> 1 191.7016 10 

> 2 178.3995 10 3 

> 3 167.3385 10 6 

fm <- lme(resp ~ as.factor(dose) :time, neurodeg, ~time|id) 
muH <- f ixef (fm) [-1] # extract estimates 
covH <- vcov(fm) [-1 ,-1] 

The estimated slopes are p, = { — 5.099, -4.581, -3.220, -2.879, -3.520) with corresponding esti- 
mated variance-covariance matrix with compound symmetry structure with diagonal elements 0.149 
and off-diagonal elements 0.0094. 

The optimal contrasts corresponding to the candidate models are calculated using the formula 
in ([3]), with S given by the estimated variance-covariance matrix of p,. The DoseFinding package 
includes the function optContr to calculate optimal contrasts based on ([3]) as follows. 
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doses <- c(0, 1, 3, 10, 30) 

mod <- Mods(emax = 1.11, quadrat ic=-0. 022, exponential = 8.867, linear = NULL, 

doses = doses) # definition of candidate shapes 
contMat <- optContr (mod , S=covH) # calculate optimal contrasts 

The MCTtest function in the DoseFinding package implements the optimal model contrast tests 
for p, based on the multiple comparison approach described in Section 12. 2[ 

MCTtest (doses , muH, S=covH, type = "general", critV = T, contMat=contMat) 

> . . . 

> Multiple Contrast Test: 

> t-Stat adj-p 

> emax 4.561 <0.001 

> quadratic 3.680 <0.001 

> linear 2.274 0.0249 

> exponential 1.277 0.1818 
> 

> Critical value: 2.275 (alpha = 0.025, one-sided) 

The Emax, quadratic and linear model contrasts are all significant, but the exponential model failed 
to reach significance. Therefore, the significance of a dose-response signal is established and we can 
move forward to estimating the dose-response profile and the target dose. 

Two approaches can be used for model fitting in this example: the two-stage GLS non-linear dose- 
response fitting method described in Section [3l or mixed-effects modeling (linear and nonlinear) 
incorporating a parametric dose response model for the progression slope fi{x). We consider first the 
two-stage GLS method, which is implemented in the f itMod function in DoseFinding, illustrated 
in the call below for the Emax model. 

f itMod(doses , muH, S=covH, model="emax" , type = "general", bnds=c(0.1, 10)) 
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> Dose Response Model 
> 

> Model : emax 

> Fit-type: general 
> 

> Coefficients dose-response model 

> eO eMax ed50 

> -5.181 2.180 1.187 

The gAIC values (as discussed in Section 13. ip corresponding to the fits for tlie Emax, quadratic, 
and linear models are, respectively: 10.66, 11.07 and 24.22, indicating the better adequacy of the 
Emax model. Note that the DoseFinding package also includes an MCPMod function that performs 
MCTtest, model selection and model fitting in one step. 

The mixed-effects model fit approach in this case is illustrated below for the Emax model using the 
nlme function 

nlme(resp ~ bO + (eO + eM * dose/(ed50 + dose))*time, neurodeg, 
fixed = bO + eO + eM + ed50 ~ 1, random = bO + eO ~ 1 I id, 
start = c(200, -4.6, 1.6, 3.2)) 

Fixed: bO + eO + eM + ed50 ~ 1 

bO eO eM ed50 

200.451303 -5.178739 2.181037 1.198791 

The estimated fixed effects from the NLME model are quite close to the estimates obtained via the 
GLS two-stage method. The AIC values corresponding to the Emax, quadratic and linear models 
under the mixed-effects model fit are, respectively: 8352.60, 8353.10 and 8365.79 confirming the 
Emax as the best fitting model. It is intriguing to see how similar the differences in AIC between 
the different models are to the differences in gAIC values. 
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Estimates for the target dose, that is, the smallest dose producing an effect greater than or equal to 
the target value of 1.4, can be obtained with either of the model fitting approaches. The resulting 
estimated target doses are 2.13 under the two-stage GLS method and 2.15 under the NLME model. 
Alternatively, model averaging methods could have been used to estimate the target dose and the 
dose-response profile. 

4.2 Binary and placebo-adjusted data 

In this section we will go through two concrete examples on how to use the DoseFinding R pack- 
age to apply MCPMod to binary data and placebo-adjusted normal data. Only the required R 
commands are given here, but not the output. 

4.2.1 Binary endpoint 

This example is based on trial NCT00712725 from clinicaltrials.gov. This was a randomized, 
placebo-controlled dose response trial for the treatment of acute migraine, with a total of 7 active 
doses ranging between 2.5mg and 200mg and placebo. The primary endpoint was "being pain free 
at 2 hours postdose", i.e., a binary endpoint. The analysis presented here is a post-hoc analysis. 

As a reasonable set of candidate models and contrasts, we select 4 different shapes of the sigmoid 
Emax model f{x,0) = Eq + EmaxX^/{x^ + ED^q), which cover a very wide range of monotonic 
shapes and a quadratic model to safeguard against the possibility of a unimodal dose response 
relationship. The Mods function is used for that, and one can also plot the candidate shapes. 

doses <- c(0,2.5,5, 10,20,50, 100,200) 

models <- Mods(sigEmax = rbind(c(2.5, 1) ,c(10, 1) ,c(50, 3) ,c(100,2)) , 
quadratic = -1/250, doses=doses) 

plot (models) 
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The first stage of tlie two-step MCPMod approacli consists of fitting a model witli ANOVA-type 
parametrization to the data to obtain estimates p, and its asymptotic covariance matrix S. The 
logistic regression model is used here, which means that the candidate models are formulated on 
the logit scale (other scales could be used). The ANOVA logistic regression model can be fitted as 
follows. 

## data from NCT00712725 study 

dosesFact <- as. factor (doses) ## treat dose as categorical variable 

N <- c(133, 32, 44, 63, 63, 65, 59, 58) 

## % of patients painfree at 2h post-dose 

RespRate <- c(13,4,5, 16, 12, 14, 14,21)/N 

## fit logistic regression (without intercept) 

logfit <- glm(RespRate~dosesFact-l , family = binomial, weights = N) 
muHat <- coef (logfit) 
S <- vcov(logfit) 

Now all subsequent inference only depends on muHat and S obtained from the logistic regression. 
The multiple contrast test from 12.21 using optimal trend contrasts can be produced as follows 

MCTtest (doses , muHat, S=S, models = models, type = "general") 

All contrasts are significant. The modeling step can now be performed using the f itMod function. 
Here for illustration we fit the sigmoid Emax model and the quadratic model. 

modSE <- f itMod (doses , muHat, S=S, model = "emax", type=" general") 

modQuad <- f itMod (doses , muHat, S=S, model = "quadratic", type = "general") 

gAIC(modSE) jgAlC (modQuad) 

A comparison of the gAIC values reveals that the sigmoid Emax model provides a better fit than 
the quadratic model. 
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Above we performed the different steps of the MCPMod procedure separately. One could alterna- 
tively have used 

MCPMod (doses, muHat, S=S, iiiodels=models , type = "general", Delta = 0.2) 
directly. 

4.2.2 Fitting on placebo-adjusted scale 

For this example we use the IBScovars data set from the DoseFinding package, taken from [3]. The 
data are part of a dose ranging trial on a compound for the treatment of irritable bowel syndrome 
with four active doses 1, 2, 3, 4 equally distributed in the dose range [0, 4] and placebo. The primary 
endpoint was a baseline adjusted abdominal pain score with larger values corresponding to a better 
treatment effect. In total 369 patients completed the study, with nearly balanced allocation across 
the doses. 

The endpoint is assumed to be normally distributed and it is of interest to adjust for the additive 
covariate gender. While the DoseFinding package can deal with this situation exactly, here we 
illustrate using the placebo-adjusted (effect) estimates. Note that, in the case of time-to-event data, 
one would proceed similarly. Here we only illustrate fitting an emax model, using the MCTtest and 
MCPMod functions is analogous to the calls in Section 14.2.11 but using the additional argument 
placAdj = TRUE. We plot the fitted model together with confidence intervals for the model fit and 
the ANOVA type effect estimates. 

data (IBScovars) 

anovaMod <- lm(resp~f actor(dose)+gender , data=IBScovars) 

drFit <- coef (anovaMod) [2 : 5] # placebo adjusted (=effect) estimates at doses 

S <- vcov (anovaMod) [2:5,2:5] # estimated covariance 

dose <- sort (unique (IBScovars$dose)) [-1] # vector of active doses 
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## now fit an emax model to these estimates 

gfit <- f itModCdose, drFit, S=S, model = "emax", placAdj = TRUE, 

type = "general", bnds = c(0.01, 2)) 
plot (gfit, CI = TRUE, plotData = "meansCl") 

5 Conclusions 

The extended MCPMod methodology, together with its corresponding software implementation in 
the DoseFinding package in R, greatly broaden the scope of application of the original MCPMod 
approach. Most type of endpoints, and associated model-based analyses, utilized in dose finding 
studies can be handled in the context of the extended approach. 

Further extensions of the approach discussed here are possible and of interest in practice. An 
increasing number of indications and drugs require regimen selection, in addition to the more 
traditional dose selection. Different approaches can be considered in the context of MCPMod, or 
its extension, discussed here. One could focus on estimating the exposure-response relationship, 
for example, combing dose and regimen into one model covariate. The much larger number of 
exposure values, compared to dose levels, could pose a problem for the derivation of optimal model 
contrasts and for the MCP step, more generally. Dose-time response modeling in the context of 
model uncertainty provide another venue for extending MCPMod. Further research is needed in 
those topics. 

Model-based dose finding methods, such as the extended MCPMod, provide better understanding 
of the dose-response relationship, generally translating into more accurate dose selection for con- 
firmatory trials. Realizing the full potential that these methods have to offer, however, requires 
changes in the way Phase II studies are traditionally designed. By and large, dose finding studies 
are planned as mini Phase III trials, using hypothesis tests to select the dose, or doses, to bring 
forward to the confirmatory stage. Relatively few doses (typically two active treatment arms and 
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placebo) are used in such Phase II studies, making it hard to entertain any type of modeling. The 
sample size derivation in this type of studies is based on power calculations for detecting at least one 
dose significantly different from placebo. The resulting number of subjects is typically inadequate 
for proper estimation of target doses (and dose response modeling). A discussion of a different 
balance in resource allocation between Phases II and III, taking into account the overall probability 
of program success, is long overdue. Utilization of larger number of doses (e.g., 4 or 5), coupled 
with larger sample sizes, would go along way in enabling model-based methods to improve dose 
selection in Phase II and, as a result, the probability of success in Phase III. 



A Derivation of Optimal Contrasts 

In this section the closed form solution for the optimal contrasts for the case of a general covariance 
structure is derived. Optimality here refers to maximum power of the univariate contrast test, if a 
specified mean vector fi (with corresponding positive definite covariance matrix S) is true, which 
means that the non-centrality parameter 

^(c,m) 



VdSc 

needs to be maximized with respect to c, subject to c'l = 0. 



Writing Cq = y — lK-i-^K~i j , the constrained maximization problem is equivalent to the uncon- 
strained maximization of -fe^;^^- This, however, is the solution to the generalized eigenvalue 
problem 

CoHf^'C'oX = ACoSCo'a;, 

see e.g. [1], formula (2.66). As Cq^^i'Cq is of rank 1, it has only a single non-zero eigenvalue. 
Thus, it is immediately clear that c = const ■ (CoSCo)"^Co/^, const 7^ is the only solution to the 
generalized eigenvalue problem. We further note that c = const - (CoSGq)~^Co/x = const - ^"^(^i — 
YT^zrrjl). It is clear that the optimal solution is invariant with respect to addition or multiplication 
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of any scalar constant to the vector /x, which is why one can also use the standardized mean vectors 
instead of n, which then gives the result in formula ([3]). 
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B Placebo- adjusted dose response modeling 



In a few cases one would like to perform MCPMod on placebo-adjusted estimates, for example when 
there are (additive) covariates in the model, or when using a Cox PH model (where one can only 
obtain control-adjusted estimates). In what follows we will first demonstrate the equality of test 
statistics, and calculate optimal contrasts. Then 



B.l Test statistics and optimal contrasts 

If we want to test a linear contrast of the responses per dose group, it does not matter whether we 
fit a placebo adjusted curve or include the placebo group as a response and then test contrasts to 
placebo. 

To see this, consider the ANOVA estimate 

n-N{tx,S) (9) 

where the first component /iq of the vector p, corresponds to the placebo response. The test 
contrasts can then be produced by multiplication of p, with the {K — 1) x K contrast matrix 
Co = ^— Ia'-i^Ia'-i^ J where Ix-i is a column vector of ones of size K — 1 and Ir-i the K — 1 
dimensional identity matrix. We obtain the corresponding contrast 

Pc-Nific^Sc) (10) 
with p,Q = Cop,, fifj = CqI-i and Sc = CqSCq. 

The contrast test statistic in model is of the form 

c^p 

mc = max —= subject to c'l^ = 0, (11) 
Vc'Sc 



with c such that nic attains a maximum. In model ffTOj) . the restriction on c is already absorbed 
in the matrix Co and the test statistic takes the form 

d'Co/x , . 

mp = max — =, (12) 

v/d'CoSC[,d 
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where d is no longer a contrast and again d is chosen such that mp is at its maximum. Now, it 
can be seen that c = Cpd: Setting c = Cgd imphes c'l^ = 0, as CqIk = 0. Hence, rric > mp. 
On the other hand, if c'Ik = 0, then there must be some d G R^~^ such that c can be written as 
c = Cgd, since the rows of the {K — 1) x i^-matrix Co provide a base of the {K — l)-dimensional 
hyperspace orthogonal to Ik in R^. Consequently, mc mp. It follows that mc = mp and that 
if d maximizes f|T2|) . then Cgd maximizes f fTTj) . 

Specifically the optimal d"^* can be calculated as d"^* = S'^^fiQ, as the sum to constraint is 
removed. 

B.2 Dose Response Model Fitting 

In the two-stage generalized least squares fitting procedure one minimizes the criterion 

{f{x,9)-nyS~\f{x,9)-n). (13) 

When we only have p,^ one would not fit a full dose response model + (^if'^i^, 0^) but work with 
a model without the intercept ^o- fcix,0) = 9if(x,6^). The optimization criterion proposed for 
placebo-adjusted is then 

{fc{x, 0) - pc)'Sc\fc{x, 9) - Tic)- (14) 

When /0(0, 6*°) = one can see that f|T3l) and f|T^ are equal. This follows from the fact that f|T^ 
is equal to 

(Co(/(a;, 0) - n)nCoSC',)-\Co{f{x, 0) - /x)), 
and because C'q{CqSC'q)^^Cq = (which follows from multiplication from the left with CqS). 
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C Proof of the result in section 



Let ^0 denote the true value of the parameter and /Xq = /(x, 6q) where aj is a known vector of 
fixed values. Assume that the following conditions are satisfied: 

(Al) There exists a symmetric, positive definite estimate S of the covariance matrix of /x with 

p 

ttnS — !• S for a positive nondecreasing sequence a„ converging to oo as n — oo, and a positive 
definite, symmetric matrix S G M'^^'^. 

(A2) If A^(0, S) denotes the multivariate normal distribution with mean and covariance matrix 
S, then an^/'^ifi — /Xq) -4 A^(0,S), where -4 denotes convergence in distribution. As a 
consequence of a„ — )• oo, the estimate /x is consistent, ^.e., /i — )■ /Xq (see e.g. Serfling, 1980, 
p.26). 

(A3) The mapping i— M'^', i— )■ f{x, 0), with a; G M'^ is a bijective function of which is twice 
differentiable in an open region around ^o- 

Under these assumptions is a consistent estimator of 0, i.e., 0^0 and y/(hi{0 — 0o) — > 
N{O,B{0oyM{0o)B{0o)), where M(6») = F{0yA'EAF{0) and B{0) = {F{0yAF{0))-\ 

Proof: 

Let ^(0) = in- fix, 0))'A„(m - fix, 0)) and *(0) = (/x - /(a;, 0))'A(^ - /(a;, 0)). 

First we note that consistency of the estimator is easy to establish using standard theory for the 
consistency of M-estimators. For example the three conditions in Theorem 5.7 from [3S] are easy 
to verify (Al)-(A3). 

The proof of the distribution of ^^(0 - 0o) A N{0, B{0Qy M{0q)B{0q)) works along the lines of 

[33| ch. 12.2.3], which we restate here with the modifications needed for our situation. 

As minimizes ^'(0), we have '^^i^^ = 0. Thus, by the mean value theorem, there is a between 
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6q and 6 such that ^ 

Hence _ 

- 0o) = (^^). (15) 

We now show that (i) a/o^^^^ is asymptotically normal, and that (ii) ^ '^dfde''' ) converges in 
probability to a non-singular matrix. 



(i) 



-2F(0)A„(/2-/(a^,0)), 



de 

where F{0) is the d x k matrix of partial derivatives ^^^^f^- {i = 1, k, h = 1, d). 
Since ya;;(Ai-/(a;,6/o)) 4iV(0,S), 

^d*(0o) ^ _2y^i.(0^) A„ O - f(x, Oo)) A iV(0, 4F{eo)A^A'F{eoy) 



de 



(ii) Differentiate the second term twice to get 

d^^(0) 



dOdO' 



-2{u - F{e)A^F{e)') 



where U is the dx d matrix with h-i\i column given by 

'11^ A (a- fix 6)) 

Now /X — ^ ^o) si'iid ^ ^ so all entries of U converge to 0. In total we get 

Defining M{e) = F{e)A'SA' F{ey and 3(6) = {F{e)AF{eyy^ and inserting the results from (i) 

and (ii) into ( ITSjl . one obtains that the asymptotic distribution of y/a^{6—6o) is A^(0, B{6oy M{6o)B{6o)). 

□ 
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Figure 6: Time- To- Event endpoint 
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Figure 8: Normal endpoint 
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Figure 9: Root mean squared error for estimating the dose response (averaged over the doses 
available), count endpoint. 
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Figure 10: Root mean squared error for estimating the dose response (averaged over the doses 
available), time-to-event endpoint. 
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Figure 11: Root mean squared error for estimating the dose response (averaged over the doses 
available), binary endpoint. 
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